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ABSTRACT 

Perturbation  theory  in  celestial  mechanics  has  a  long,  rich  history  of 
failure  stretching  back  to  Newton.  I  believe  that  the  causes  of  this  are 
two- fold.  One  problem  is  the  difficulty  of  dealing  with  the  mathematical 
structure  used  in  celestial  mechanics  to  express  perturbation  theory  as 
opposed  to  the  constructs  used  in  field  theories  (eg.  trajectory  equations  vs. 
linear  second  order  partial  differential  equations).  The  second  flows 
directly  from  this  and  relates  to  the  misapplications  of  certain  mathematical 
techniques  (averaging,  series  expansions)  within  the  context  of  perturbation 
theory.  These  incorrect  analyses  usually  appear  in  second-order  theories  such 
as  Kozai's  (1959)  artificial  satellite  theory.  Ideally  this  report  would 
clearly  illustrate  the  nature  of  these  difficulties  utilizing  a  complex  (but 
exactly  soluble)  physical  model  intimately  tied  to  the  two-body  problem  and 
then  go  on  to  lay  the  foundations  for  a  new  perturbation  theory.  I  believe 
that  that's  exactly  what  is  accomplished  herein  except  that  the  hints  of  the 
base  of  this  new  mathematical  formalism  are  severely  limited.  The  exactly 
soluble  physical  model  is  the  three  dimensional  harmonic  oscillator  compli¬ 
cated  by  anisotropic  terms,  anharmonic  terms,  and  air  resistance.  The  deep 
connection  is  provided  by  Bertrand's  theorem  which  is  also  proved. 
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I.  INTRODUCTION 

There  is  no  problem  In  celestial  mechanics  that  perturbation  theory 
provides  a  satisfactory  solution  to  (of  which  I  am  aware).  Saturn's  rings, 
the  Kirkwood  gaps,  Jupiter's  satellites,  the  Sun-Jupl ter- Saturn  system,  etc., 
all  exhibit  resonances,  unusual  structures,  or  some  other  peculiarity  not  yet 
deduced  from  celestial  mechanics.  Celestial  mechanics  may  provide  explanations 
(as  in  the  fact  that  the  Kirkwood  gaps  occur  at  simple  multiples  of  Jupiter's 
mean  motion)  but  not  predictions.  Why  Is  this  state  of  affairs  so  poor?  I 
believe  that  the  reason  is  two-fold  (at  least).  First  of  all  celestial 
mechanics  is  a  trajectory  theory  and  not  a  field  theory.  The  field  theories 
with  good  perturbation  techniques  rely  on  linear  second  order  differential 
equations  not  non-linear  first  order  ones  (eg.  the  Hami 1  ton- Jacobi  equation). 
Secondly  the  mathematical  nature  of  the  problem  in  celestial  mechanics  is 
compounded  by  a  poor  choice  of  coordinate  systems.  It  does  no  good  to  express 
perturbation  theory  in  terms  of  the  orbital  element  set,  because  we  will  then 
immediately  and  intuitively  understand  the  result,  if  we  can  never  solve  any 
problems  in  this  coordinate  system.  This  becomes  compounded  when  our  frustra¬ 
tion  with  our  failures  further  leads  us  astray  in  the  use  of  unjustifiable 
analytic  techniques.  Examples  of  the  latter  include  incorrect  power  series 
expansions  and  the  method  of  averaging. 

It  is  possible  that  the  reader  Is  not  as  aware  of  the  shortcomings  and 
pitfalls  of  the  misapplications  of  perturbation  theory  in  celestial  mechanics 
as  the  above  paragraph  suggests  Is  the  case.  This  report  will  serve  to  inform 
the  reader  of  the  situation  by  a  carefully  constructed  sequence  of  illustra¬ 
tions.  As  the  nature  of  the  demonstration  is  by  counter-example,  it  follows 


that  each  counter-example  must  in  itself  be  exactly  soluble.  Moreover  it  is 
clear  that  no  such  meaningful  example  is  likely  to  come  from  the  normal  realm 
of  problems  dealt  with  in  celestial  mechanics.  However  it  would  be  preferable 
to  deal  with  a  system  as  close  to  the  two-body  problem  as  possible.  The 
existence  of  Bertrand's  (1873)  theorem  is  used  to  find  the  suitable  connection 
and  the  alternative  physical  model--the  three  dimensional,  isotropic,  simple 
harmonic  oscillator. 

Bertrand's  theorem  simply  states  that  of  all  possible  central  force 
potential  functions  only  two  give  rise  to  bounded,  closed  orbits.  One  is 
the  two-body  potential  (- k/r )  and  the  other  is  the  isotropic,  simple  harmonic 

o 

oscillator  (kr  /2).  The  next  section  follows  Arnold's  (1978)  demonstration 
of  this  result.  Following  that,  motion  under  the  three  dimensional,  isotropic 
simple,  harmonic  oscillator  potential  is  discussed  (the  orbits  are  ellipses 
with  force  center  at  the  center).  A  perturbation  theory  for  these  orbits  is 
developed  and  a  set  of  first  order  differential  equations  is  derived.  These 
are  exactly  analogous  to  Lagrange's  variation  of  parameter  equations  in  the 
two-body  problem. 

The  fourth  Section  deals  with  the  mathematical  foundations  of  perturba¬ 
tion  and  an  attempt  to  use  the  variation  of  parameters  equations  when  a 
perturbation  yielding  anisotropy  is  applied.  This  attempt  very  closely 
follows  Kozai's  (1959)  artificial  satellite  theory  in  format.  The  purpose 
is  not  to  criticize  Kozai  but  rather  to  illustrate  the  invalidity  of  a 
technique  widely  used  in  celestial  mechanics.  This  section  concludes  with 
the  development  of  a  new  set  of  variation  of  parameters  equations  utilizing 
a  different  set  of  constants  (that  is  other  than  the  orbital  element  set). 


In  this  representation  the  perturbation  equations  are  exactly  soluble  to 
all  orders— indeed  the  whole  system  is  exactly  soluble  in  closed  form.  The 
misleading  inferences  of  the  classical  approach  can  be  clearly  seen  now. 

The  fifth  section  introduces  anharmonic  perturbations.  These  too,  in 
the  right  basis,  present  a  set  of  perturbation  equations  for  which  we  can 
obtain  the  solution  to  all  orders  (explicitly)  as  well  as  in  closed  form. 

This  further  elucidates  the  nature  of  the  failings  of  other  approaches. 

In  the  final  section  I  discuss  non-conservative  perturbations  due  to 
drag  and  "third-body"  forces  for  the  harmonic  oscillator.  This  is  very 
short  as  are  some  remarks  concerning  additional  representations.  By  now  I've 
proved  my  points  and  our  energies  should  be  devoted  to  productive  rearrange¬ 
ments  of  perturbation  theory  for  the  two-body  problem. 


II.  BERTRAND'S  THEOREM 

According  to  Plummer  (1908),  in  1873  Bertrand  proved  that  the  inverse 
square  law  force  and  the  direct  linear  force  law  are  the  only  two  central 
force  laws,  derivable  from  a  potential  function,  which  yield  bounded  closed 
orbits.  On  pages  2-9  of  his  book  Plummer  presents  a  "proof"  of  this  result. 
Goldstein's  second  edition  of  Classical  Mechanics  (1982)  reproduces  Plummer's 
result  in  modern  language.  It  does  not  appear  to  me  that  this  proof  is 
rigorous.  Arnold  (1978),  without  ever  mentioning  Bertrand,  does  provide  the 
outlines  of  a  rigorous  proof  of  this  result.  While  I  shan't  fill  in  all  of 
the  analytical  steps  for  the  reader,  immediately  below  I  flesh  out  Arnold's 
logic.  The  importance  of  the  result  for  this  discussion  is  that  it  yields 
a  deep  connection  between  the  two-body  problem  and  the  three  dimensional, 
isotropic,  simple  harmonic  oscillator  problem.  (I  felt  it  incumbent  upon  me 
to  debunk  the  misapplications  of  perturbation  theory  via  a  physical  situation 
either  identical  to  the  1/r  potential  or  intimately  associated  with  it.) 

The  logic  of  Arnold's  argument  can  be  simply  outlined.  First  of  all, 
if  all  bounded  orbits  are  closed,  then  since  a  circular  orbit  is  always 
possible  for  a  central  force,  all  nearly  circular  bounded  orbits  should  be 
closed  too.  Moreover  this  (eg.  being  closed)  must  be  true  independent  of  the 
radius  of  the  circular  orbit.  The  implications  of  these  statements  are  that 
the  potential  energy  is  either  a  power  law  U(r)  =  Arp  (pj>-2,p/0)  or  a 
logarithmic  function  U(r)  =  B Znr  (the  p=0  case).  By  considering  the  con¬ 
straints  of  "closedness"  on  the  apsidal  angle  Arnold  then  shows  that  the 
logarithmic  case  can  not  be  closed  while  there  are  two  possibilities  for 
the  power  law  case.  These  two  deal  with  positive  power  laws  (p>0)  and  the 
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meaning  of  bounded  and  with  negative  power  laws  which  yield  infinitely  tight 
binding.  After  an  examination  of  these  two  limiting  cases  it  turns  out  that 
only  p=+2  and  p=-l  are  left.  QED. 

The  conservation  of  energy  equation  for  a  central  force  after  the  usual 
reduction  to  the  equivalent  one-body  problem  (cf.  Eq.  2  below)  yields 


2  2 

Hence  the  "pseudo-potential  energy"  V  *  U  +  L  /2mr  . 


From  this  we  have  that 


Jdt  =  +  Jdr  /  [2(E-V)/m]1/2 
2* 

and  from  L  =  mr  0, 


(l/r2)dr 

[2m(E-V)]1/2 


Orbits  exist  such  that  E  _>  V(r)  for  a  fixed  value  of  L.  This  inequality 
gives  one  (or  several  depending  upon  V)  regions  with 


0  <  r  .  <  r  <  r  <  « 

-  min  -  -  max  - 


Clearly  if  0  <  r  .  <  r  <  °°  then  the  motion  is  bounded  (r  w  <  »)  and 

—  min  max  max 

also  takes  place  in  the  annulus  with  inner  radius  r  .  and  outer  radius  r  . 

mi n  max 

As  9  increases  r  varies  between  these  two  limits.  The  angular  extent  between 
successive  apsides  (extrema  of  r)  is  called  the  apsldal  angle.  It  is  given 
by  the  expression 


(1) 


A  closed  orbit  is  one  such  that  0  is  commensurate  with  2tt.  Finally  in  the 


special  case  when  E  is  equal  to  a  minimum  of  V  then  rm.  =  rmax  =  R  and  the 
orbit  is  a  circle  of  radius  R.  Now  that  the  stage  has  been  set, the  argument 
can  be  filled  in. 

First  let  us  calculate  the  value  of  0  for  a  nearly  circular  orbit.  Since 
the  orbit  is  nearly  circular  E  must  be  slightly  larger  than  the  minimum  value 
of  V.  Near  r=R  V  can  be  written  as  (since  it's  a  minimum  V ' (R)  =  0,  V"(R)>0) 

2 

V (r )  =  V (R )  +  V"(R) 

The  extremes  of  r  are  found  by  solving  E-V(r)  =  0.  They  are 
R  +  {2[E-V(R)]/V"(R)}1/2 

with  rmax  given  by  the  plus  sign  and  rmin  given  by  the  minus  sign.  Use  the 
definition  in  Eq.  (1)  of  0,  plus  the  change  of  variables 


so 


0riw  -  TrL/R2/V"(R)m 


In  terms  of  U, 


0circ  “  ir[U7(3U'+RU")]1/2/vfir 

Now  impose  the  fact  that  we  want  this  result  to  be  independent  of  the 
value  of  R.  This  simply  means  that 


U' 

3U'+RU-1 


=  constant 


Choose  the  form  of  the  constant  to  be  l/(2+p)  and  integrate.  One  gets 
directly  that 

U  =  Ar^  (p>-2,pj*0)  or  U  *  Btnr 


within  an  arbitrary,  unimportant,  additive  constant.  From  the  formula  for 
0circ  we  now  get  that  ^set 

©circ  3  ir/ (2+p)1  (p=0  is  the  logarithmic  case) 


This  completes  the  first  part  of  the  result. 

If  p=0  (the  logarithmic  case)  then  0circ  =  ir/,/2  which  is  not  a  rational 
multiple  of  2m.  Hence  this  functional  form  can  be  dismissed  from  further 
consideration.  If  p>0  then  as  r-**>,  U(r)-*».  Therefore,  to  remain  bounded,  E-*» 
too.  The  general  limit  as  E-*»  of  0  Is  tt/2 .  To  see  this  make  the  change  of 
u=L/r  in  Eq.  (1).  Then 
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where  W(x)  =  U(L/x)  +  x^/2m.  Further  let  y  =  x/umav  so 

niaX 

1 

0  =  /  dy/{2m[w(l )-w(y)]}1^2 
y  . 

where  w(y)  =  y2/2m  +  U(L/y  %,ax)/u2max.  Clearly  as  E-*-  »,  umax->  »,  ymin  -►  0 

and  the  last  term  in  w  is  negligible.  The  tt/2  result  follows  after  a  trivial 

integration.  Now  if  tt/(2+p)^2=tt/2  for  p>0,  p=2  and  we  have  the  isotropic 

harmonic  oscillator  as  a  possibility. 

The  last  piece  involves  p<0.  Here  there's  the  possibility  that 

as  r+O.  First  we  show  that  in  general  as  E-*-°°,  0=ir/(2-q)  for  power  law 

potential  energies  U  =  -kr”q,  qe(0,2).  The  proof  is  as  above  and  leads  to 
1 

0  =  f  dx/(xq-x2)1/2  =  Tr/(2-q) 

0 

(This  is  one  of  several  misprints  in  this  section  of  Arnold's  book.)  Now 
replace  -q  by  p  and  demand  equivalence  with  the  circular  orbit  result  or  that 

ir/(2+p)  =  Tr/(2+p)1^2 

Only  p  =  -1  satisfies  this  case--the  inverse  square  law. 


III.  THE  ISOTROPIC,  THREE  DIMENSIONAL,  SIMPLE  HARMONIC  OSCILLATOR 

Just  to  be  sure  that  we  understand  the  terminology  the  potential  energy 

2 

is  U  =  kr  /2.  The  force  constant  is  k  and  it  is  the  same  in  all  directions 
(isotropic).  If  the  force  constant  differed  in  different  directions,  then 
U  =  r*fe'r/2  where  fe  is  now  a  diagonal  tensor 


This  would  be  an  anisotropic,  three  dimensional,  simple  harmonic  oscillator, 
p  o  2  2 

Since  r  =x  +y  +z  this  explains  the  three  dimensional  aspect.  The  simple 
harmonic  part  is  due  to  the  fact  that  U  is  quadratic  in  x,y,z.  The  resulting 
motion  is  then  sinusoidal  with  a  period  independent  of  the  amplitude.  The 
next  Section  discusses  anisotropic  oscillators  and  the  one  following  that 
anharmonic  ones. 

The  force  is  £  =  -  VU  *  -kr  =  mr.  Energy  conservation  is  obvious  and 
we  set 

E  =  U+  |r*r>0  (2) 

Furthermore  since  £  is  central  L  =  rxmr  isa  constant  vector.  Therefore 
the  orbit  lies  in  a  plane  given  by  r  '  £  =  0  and  the  problem  can  be  reduced 
to  a  two  dimensional  one.  Set 


£  =  L  (sini  sinfl,  -  sini  cosn,  cosi) 


(3) 


where  i  (the  inclination)  and  ft  (the  longitude  of  the  ascending  node)  specify 
the  orientation  of  and  L  is  its  magnitude.  Finally  introduce  polar  coordi¬ 
nates  r,0  in  the  orbital  plane  so  that  the  equations  of  motion  are 


m(r  -  r02)  =  f(r)  =  -dU/dr  =  -kr 

?  Jt  <*>  ■  0 

A.  The  Orbit 

As  usual  introduce  Binet's  transformation  u  =  1/r  and  switch  from  time 
derivatives  to  derivatives  with  respect  to  0  via  (prime  =  d/d0) 

A 

r  =  r'0  =  -u'0/u c  =  -Lu'/m 


Then  we  find  that  the  equations  of  motion  are 


u"  +  u 


(4) 


Define  the  mean  motion  n  by 
n2  =  k/m 

A  single  integration  of  Eq.  (4)  gives 


(5) 


Comparison  with  the  conservation  of  energy  equation  (2)  in  the  u  form  shows 

2 

that  the  constant  is  2Em/L  .  Hence, 


and 


o  5  _ L  /m 

rmax  =  a  E-(E2-L2n2)] 


E+(E2-L2n2)^2 


.2  _  l2  _  .2/-,  .2x  _  LVm 

r^"  ‘ b  ' a  °‘e  ’ ' 


.meW? 


At  r  a  a  or  b  r  =  0  and  r‘  *  0  hence  their  Identification  with  the  extrema  of 
r. 

Suppose  that  the  motion  starts  at  0  =  0Q  where  r  *  a.  Then  thereafter 
r'  <  0  so 


0 

+  J  d<fc 


Without  loss  of  generality  choose  0QSO  and  let  another  orientation  angle,  w, 
define  the  line  of  the  major  axis  ("the  argument  of  peri  apse").  Then 


is  the  equation  of  the  orbit.  If  E  s  rcos©,  n  s  rsin©  are  rectangular 


coordinates  in  the  orbital  plane,  then  this  is  equivalent  to 
S2/a2  +  n2/b2  =  1 

This  is  the  equation  of  an  ellipse,  center  at  the  origin,  major  axis  a,  minor 
axis  b  (and  therefore  eccentrlcy  e),  whose  major  axis  (the  line  of  apsides) 
lies  on  n  =  0. 

B.  The  Time  Dependence 

Return  to  the  conservation  of  energy  equation  (2)  and  solve  for  r2, 

r2  =  \  (a2-r2)(r2-b2) 
r 

Suppose  r=a,  0=0  corresponds  to  t=T  ("the  time  of  periapse  passage").  Since 
a>  r,  r<  0  for  t>T  so 


/ 


[(a2-x2)(x2-b2)]1 


or 

r2  =  a2cos2M  +  b2  s1n2M  (8) 

where  the  mean  anomaly  M=n(t-T).  Note  that  9  Is  the  true  anomaly  and  the 

"eccentric  anomaly"  too. 

2  2  2 

However  r  =  £+rf  so 


C  =  acosM,  n  *  bsinM 


(9) 


Note  too  that  L  =  mr  8  *  2mdA/dt  where  A  Is  the  areal  rate.  Whence, 


L  =  2m(irab/P)  =  mnab  (10) 

where  the  period  is  P  =  27r/n. 

To  determine  the  time  dependence  of  0  use  L=mr29  and  the  expression  (8) 
for  r2(t),  viz. 


_ dr _ 

"5  5  5  7~ 

a  cos  y+b  sin  y 


,  y=n(x-T) 


The  result  is 

tane  =  (1-e2)^2  tanM,  M=n(t-T) 

This  is  "Kepler's  equation". 

To  summarize  the  three  dimensional  motion  is  given  by 

(sin9\  /  cosM  \ 

cos0  j  55  aSm-e2)1^2sinMj 


(ID 


(12) 


where  the  unitary  rotation  matrix  S^R^-njR-i (-i )R3(-u).  The  R  matrices  are  the 
usual  elementary  rotation  matrices.  The  angular  momentum  has  already  been 
given;  the  energy  E  is 


C.  A  Perturbation  Theory 


Now,  because  we  want  to  consider  more  complicated  physics  but  our 
analytical  capabilities  aren't  up  to  solving  the  more  realistic  problems, 
we  need  a  perturbation  theory.  Since  our  knowledge  of  analytical  geometry  is 
almost  intuitive  we'll  feel  more  comfortable  if  the  perturbation  theory  is 
expressed  in  terms  of  the  orbital  element  set  &  =  (a,e,w,i ,fi,T)  or  some 
function  thereof.  Hence,  confining  ourselves  to  a  conservative  disturbing 
force  -HrflR  for  the  moment,  we  derive  the  perturbation  equations  equivalent  to 

mr  =-  mVU  +  mVR 

where  U  is  the  zero'th  order  potential.  So  set  r  *  r(a^t)  and  derive  that 
3r 

7,  af  •  i-  7R 

Along  the  way  I've  Imposed  the  condition  of  osculation  on  r, 
r  =  3r/3t  or  Va  r  *  a  =  0 

w  a  —  — 

so  that  both  locations  and  velocities  can  be  computed  from  the  usual 
(eg.  unperturbed)  formulas.  Finally  we  rearrange  all  of  this  via  Lagrange 
brackets  to  read 


[arak^ak  = 


aR/Sa^ 


(14) 


where  the  Lagrange  bracket  is  defined  by 


tvV  *  (ft  %  -  fit  Hi)  (15> 

We  note  that  dCa^.a^l/dt  =  0  so  that  the  Lagrange  brackets  can  be 

evaluated  at  any  convenient  place  In  the  orbit.  I  choose  M=0  (t=T).  We 

further  note  that  there  are  (|)=  15  Independent  brackets  because  a  Lagrange 

2 

bracket  Is  antisymmetric  and  the  most  there  could  be  Is  6  =36.  The  last 
point  of  Interest  Is  to  remember  Eq.  (12)  for  r  and  realize  that  there  are 
three  types  of  Lagrange  brackets: 

Type  I:  both  a&,ak  are  one  of  1,w,n  of  which  there  are 
(2)=  3  such. 

Type  II:  a^  Is  one  of  1*oj,ft  but  a^  Is  one  of  a,e,T  (or  M)  of 
which  there  are  3*3  *  9  such  and 

Type  III:  both  a^.a^  are  one  of  a,e,T  (or  M)  of  which  there 
are  (2)  ’  3. 

So  far  this  Is  (formally)  just  like  the  two-body  problem.  In  fact  all 
of  the  brackets  not  Involving  a  are  the  same  as  therein  and  the  whole  set  of 
non-zero  results  Is 

[a.M]  *  -an(2-e2)  *  -[M,a] 

[a, ft]  =  -2bncos1  *  -[ft, a] 

[a,u>]  3  -2bn=  -[u.a] 

[e,M]  =  a2ne  -  -[M,e] 


[e.ft] 


2 

a  necosl 


(l-e‘) 


Cn.e] 


Ce  ,cj3 


(l-e*)1 


-C«.e] 


[1,ft]  *  abnslnl  *  -[ft,l] 


Alternatively  [T,e]  =  a2n2e,  [T,a]  s  an2e2 


The  last  step  is  to  put  these  results  into  Eq.  (14)  and  explicitly 
solve  for  a.  One  finds 


da  1  3R  (l-e2)1/2  3R 

dt  "  ~  7  aff  '  - * -  3u> 

ane  ,__2 
ane 


a  ne 


7  T 
a  ne 


do) 


(16) 


dt  " 

di  _ 
dt 

dO  _ 
3t  ' 


(l-e2>1/2 

^ y 

|(2-e2)  || 

+  ae 

3R 

3a  ' 

csci 

^cosi 

3R 

na2(l-e2)^2 

3w 

. 

3ft/ 

csci 

3R 

iTa'V?)1^ 

3T 

1-e 


ai 


dM  _  -1  3R  2(l-e2)  3R 

dt - ?  3a  "  7  y  3e 

nae  a  ne 


s 

I 

1 

t) 


e 


1 


a 


3 
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IV.  THE  ANISOTROPIC,  THREE  DIMENSIONAL,  SIMPLE  HARMONIC  OSCILLATOR 


Before  delving  into  the  main  subject  of  this  Section  it  seems  prudent  to 
me  to  review  the  mathematical  basis  of  perturbation  theory  within  celestial 
mechanics  and  its  relatives  in  field  theory.  This  will  clearly  reveal  the 
difficulties  which  those  of  us  who  work  in  celestial  mechanics  face  while 
demonstrating  our  collective  blindness  to  the  "right"  coordinate  system  (only 
Vinti  [1961] has  overcome  this  but  in  a  trivial  way).  Then  the  problem  will 
be  completely  solved,  in  a  variety  of  ways,  in  the  right  coordinate  system. 

A.  The  Mathematical  Foundations  of  Perturbation  Theory 

Within  classical  mechanics  we  deal  with  equations  of  motion  of  the  form 

F  =  ma 

If  we  set  v  =  dir/dt,  then  an  equivalent  first  order  system  is 

v  =  dr/dt,  £  =  mdv/dt 

If  F  is  conservative,  F  =  -W,  then  we  have 

v  =  dr/dt  ,  -VV  =  mdv/dt  (17) 

Whenever  V  is  such  that  we  can't  solve  the  problem  analytically,  we  try  to 
separate  it  into  two  parts.  One  part,  represented  by  U,  is  such  that  Eqs. 
(17)  are  exactly  soluble,  while  the  other  part  (represented  by  -R  in  the 
standard  sign  convention)  makes  the  problem  intractable.  The  solution  of 

v  =  dr/dt,  -  VU  =  mdv/dt 

has  six  arbitrary  constants  associated  with  it,  say  a.  The  parametrization 
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MmA. 


of  the  solution,  r  =  r  (a^t),  is  then  used,  as  discussed  above,  to  formulate 
Since  |R|  is  small  in  some  sense  we  can  rewrite  it  as  cR  where  e  is  a 
small  parameter.  Hence  we  are  led  to  study  the  solution  of  equations  of  the 
form  (or  systems  of  such  equations) 

dy/dx  =  €f(x,y)  ,  y(xQ)  =  yQ  (18) 

The  development  below  is  independent  of  the  dimensionality  of  y. 

I  know  of  one  rigorous  method  of  approximately  solving  Eq.  (18)— Picard's 

method  of  successive  substitutions.  Glossing  over  the  mathematical  details 

(see  Ince  1927)  we  construct  the  sequence  |yn(x)|, 

x 

yn(x)  =  V€/  fC*»yn-l^dx',  n  =  1,2,... 
x0 

The  sequence  |yn-yn_-||  converges  uniformly  and  absolutely  to  a  function  y(x) 
which  is  continuous  and  satisfies  Eq.  (18).  The  solution  of  Eq.  (18)  is 
unique  and  stable  with  respect  to  volume  perturbations  (eg.  small  changes  in 
f)  and  surface  perturbations  (eg.  small  changes  in  yQ).  That's  it  in  a  nut- 
shel 1 . 

Computationally  carrying  out  the  Picard  process  is  enormously  difficult 
in  general.  To  my  knowledge  it  has  never  been  carried  out  beyond  n=l  in 
celestial  mechanics.  The  reasons  for  this  are  two:  the  orbital  element  set 
is  an  extraordinarily  poor  reference  frame  which  exploits  none  of  the 
symmetries  of  R  (or  U)  and  the  interesting  problems  of  celestial  mechanics 
(atmospheric  drag,  oblateness  perturbations,  and  third-body  effects)  are 
not  simple.  In  an  effort  to  deal  with  these  difficulties  other  forms  of 
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approximation  have  been  tried.  One  is  to  try  and  obtain  a  solution  to  Eq. 


(18)  in  a  power  series  in  (.  For  instance  one  would  calculate  y-j  correctly 

x 

Y-j (x)  =  y-j(x)  =  yQ+ef  f(x1,yQ)dx1 

x0 

but  instead  of  computing  y2  from 

x 

y2(x)  *  y0+€/  f[x1,y1(x1)]dx1 
x0 

replace  y1  by  its  expression, 

x  X-j 

y2(x)  =  y o+c/  f£We/  f(x2»yo^dx2^dxi 
xo  xo 

and  then  expand  f  in  a  power  series  ine(9f(x,y)  *  3f(u,v)/9v|  ) 

U“ A  9 V 

X  XJ 

y2(x)  =  y0+«/  {^(x1,y0)  «af(xi  »y0)[7  f(x2*yo^dx2 j  }dxi 

xo  xo 

X  x  x1 

f2(x)  =  y0+«/  f  (x-j  ,yQ)dx-j  +  €  J  3f(x-|,yg)  £ J~  f  (x2»yg)dx2  J  dx-j 

x0  x0  x0 

If  one  continues  this  policy  to  higher  orders— expanding  f  at  each  step  so 
that  for  Y3 
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y3(x)  =  y0+f/  ^Cx-j .y2(x-, )]dx1 
x0 


x  X-j 

=  y0+£/  f{xT  »y0+€J  fCx2,yl  (x2^dx2}dx1 

x0  XQ 


=  y0+6j  (f(xry0)  +^(xry0){  /  f[x2,y1(x2)]dx2^dx1 
x0  x0 

x  x  x^  x2 

■cj"  f (x-j  »yg)dx-j+c  j  df(x-,  ,y0)  {  J  ^[x2»yQ"*"f^  ^^x3’^o^x3^x2}' 
x0  x0  x0  x0 


x  x  x1 

^X1 »yo^xl+€  J  (x-j  »yQ)  y  J  {^^x2’^0^ 
x0  X0  X0 


+  €af(x2,y0) 


2  \ 

[I  f  <x3*y0)dx3  ]  }dx2  )dxl 
x0 


a  Y0  +  €3  term 


or,  in  general. 


v*> 


pA  Ap  ^ 

■  vn-l+t  I  df(xl  ,y0*"'(j  ^V 


l"(Vy0,dxn)  dxn-l)  dxn-2’*,dxl 


•y0>(/"f<Vl 

x0 


--then  one  can  prove  that  the  sequence  (Yn>  converges  absolutely  and 
uniformly  to  a  continuous  function  of  x  that  satisfies  the  initial  condition. 
The  series  does  not  satisfy  the  differential  equation.  Hence  it  is  useless. 

If  instead  of  substituting,  expanding,  substituting  etc.  one  substitutes 
and  then  expands  one  gets  a  different  result  (beyond  y2  or  Y2).  For  instance 
return  to  y3> 

x 

y3(x)  =  y0+€J  f[x1,y2(x1)]dx1 
x0 

x  x1 

=  y(M  f{xi,yo+€J  ^Cx2*yi  (x2)]dx2}dxi 

xo  xo 

X  X-|  X2 

=  V€1  f{xl,y0+€J  f^x2,y0+€j  f(x3*y0)dx3^dx2}dxl 

x0  x0  X0 


Now  expand  from  the  inside  outwards 


i 


2 


A  A-j 

y3(x)=;  yg+€f  ^{xi  *^o+^f  [f(x2 »yo)+€af(x2,yo) 
xo  x0 


( |  f^x3’y0^dx3)]  dx2  }dxl 
x0 


x  x-|  x1  x2 

=  y0 +€I  f{xl’yO+fJ  f(x2’y0^dx2+£2J  af<x2’y(P  (|  f(x3*y0^dx3)  dx2}dxl 
x0  x0  X0  x0 

x  4  X1 

“y0+€f  pxryO)+  [€J  f(x2*y0)dx2 

xo  xo 

X1  x2 

+  £ZJ  ^f(x2,y0)  ( J  f(x3>y0)dx3^  dxg]  af(xvy0) 

x0  x0 

2  x-j  2 

+  2  f(x2»yg^dx2j  ^  ^(X-j  jYq)  |  ^x] 

X0 

through  all  terms  of  order*2  (<?2f(x,y)  =  d2f  (u,v)/«Jv2 1  ). 

u-x,v-y 

This  is  not  equal  to  y3  or  Y3  but  the  technique  represents  a  perfectly  accept 
able  approximation  scheme  (on  the  surface).  Note  too  that  the  Yn  expansion 
only  required  that  df/dy  exists  while  this  will  require  that  all  partial 
derivatives  of  f  with  respect  to  y  exist  and  be  continuous. 

The  method  of  averaging  is  yet  another  way  around  the  difficulty  which 
is  not  rigorous  either  (so  far  as  I'm  aware).  Indeed  in  the  next  subsection 
the  use  of  it  on  the  harmonic  oscillator  will  illustrate  it's  shortcomings. 
Now  what  does  perturbation  theory  for  a  scalar  field  look  like?  Suppose 


that  we  are  dealing  with 


+  (k  -U0-XU)ip  =  0 


where  the  unperturbed  problem  is 
V21J>  +  (k2-U0)ip  =  0 

Let  the  solutions  be  the  complete  set  of  orthonormal  eigenfunctions  un 

o 

corresponding  to  the  eigenvalue  kp  (degeneracy  is  an  annoyance  not  dealt 
with  here).  We  found  these  eigenfunctions  by  exploiting  the  symmetries  of 
Uq, not  by  restricting  ourselves  to  the  same  coordinate  system  for  all 
perturbations.  The  approximations  to  the  perturbed  eigenfunctions  and 
eigenvalues  are  given  by 

*(1)  ■  u  ♦  x  £  -^up 

n  p*>  *„-*p  p 

<k2>0)  '  kn  +  xurn 


♦(2>  •  ur  *  x  £  Vp  ,  +  v 

*  (kn+XUnm-kp> 


2  y  uuu 
Z-  pq  qn  c 


(k2)<2>  -  k2  +  xu  ♦  X2  £  Veil 
n  nn  ,x 


P*n  kn’kp 


and  for  the  a'th  order 


<(a)  -  u„  +  x  E  — •  • 

n  pfn  [(kz)(a_1,-k2] 


+  Xa  £ 


U-j.*  *  •  U,  u 
P5 _ znj 


w"^n  (k?-k2)(k2-k2)...(k2-k2J 


I've  listed  these  formulas  directly  from  Morse  and  Feshbach  (1953)  which  the 
interested  reader  should  consult.  In  general  these  series  converge.  Note 
that  here  the  work  consists  of  calculating  the  matrix  elements  of  U  and 
summing  the  series  whereas  in  celestial  mechanics  the  analytical  work  never 
ends. 
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B.  The  Method  of  Averaging 

Kozai  (1959)  exploits  the  method  of  averaging  in  order  to  solve  the 
main  problem  of  artificial  satellite  theory.  This  is  defined  as  to  solve 
for  the  motion  of  a  particle  under  a  potential  which  includes,  in  addition 
to  the  usual  1/r  term,  the  first  three  oblateness  terms  assuming  axial 
symmetry— namely  the  Jg.Jg  an<*  ^4  ternis-  Various  developments  of  the 
method  of  averaging  (Lorell,  Anderson,  and  Lass  1964;  Lorell  and  Liu,  1971; 
McClain  1977)  all  refer  back  to  Kryloff  and  Bogoliuboff  (1947)  and  Bugoliubov 
and  Mitropolsky  (1961)  for  a  basis.  If  this  method  has  a  rigorous  foundation 
I  can't  find  it.  As  an  example  consider  the  following  (pgs.  39-41  of  the 
last  mentioned  reference): 

We  now  go  over  to  the  discussion  of  a  method  for  developing 
the  first  instance,  asymptotic  approximations,  for  the  case  of 
oscillations  defined  by  a  differential  equation  of  the  form: 


where  e  is  a  small  positive  parameter.  We  can  arrive  at  the 
correct  formulation  of  this  method  if  we  start  from  the  physical 
concepts  defining  the  character  of  the  oscillatory  process. 

When  perturbation  is  absent  i.e.  when  e  =  0,  the  oscillations 
will,  evidently  be  purely  harmonic, 

x  =  a  cos  ip 

with  a  constant  amplitude  and  a  uniformly  rotating  phase  angle: 
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3T  =  31t  *  w 

(the  amplitude  a  and  the  phase  0  of  the  oscillations  will 
be  constants  over  time,  depending  on  the  initial  conditions) 
The  existence  of  non-linear  perturbation  (e^O)  results 
in  the  appearance  of  overtones  in  the  solution  of  equation 
(1.1),  a  factor  that  establishes  dependence  between  the 
instantaneous  frequence  and  the  amplitude,  and  finally 
gives  rise  to  a  systematic  increase  or  decrease  in  the 
amplitude  of  the  oscillations,  depending  upon  whether  the 
energy  is  expelled  or  absorbed  by  the  perturbing  forces. 

All  these  effects  are  obviously  in  the  limiting  case  (e=0). 

With  all  this  in  view  we  shall  seek  a  general  solution 
of  equation  (1.1)  in  the  form, 

2  3 

x  =  acosijj  +  eu-|(a,ij/)+  e  u2(a,i|>)  +  e  u3(a,ip)  +  ... 

Here  u-j(a,^),  U2(a,<|>)...  are  periodic  functions  of  the  angle 
tp  with  a  period  2ir  and  the  quantities  a,ip  are  functions  of 
time  defined  by  the  differential  equations: 

gl- *  GA-|(a)  +  e^Agta)  +  ...» 
if  *  go  +  eB, (a)  +  e^B«(a)  +  ... 


We  are  to  choose  suitable  expressions  for  the  functions 


u-|(a,tj;),  u2(a,^)...,  A-|(a),  B-j (a ) •  A2(a),  B2(a),...  In  such  a 
way  that  the  equation  (1.2),  after  replacing  a  and  ip,  by  the 
functions  defined  In  equations  (1 .3),  would  serve  as  a  solution 
of  equation  (1.1). 

As  soon  as  this  problem  Is  solved  and  explicit  expressions 
for  the  coefficients  of  expansions  occuring  in  the  right-hand- 
sides  of  equations  (1 .2)  and  (1.3)  are  obtained,  the  problem  of 
integrating  equation  (1.1)  is  reduced  to  that  of  integrating 
equations  (1 .3)  which  have  separable  variables,  thus  making  the 
Investigation  possible  with  the  help  of  well  known  elementary 
methods. 

We  note  that  In  the  case  represented  by  equation  (1 .1)  we 
might  establish  the  convergence  of  expansions  (1 .2),  (1 .3)  under 
very  general  conditions  for  the  function  f  • 

However,  since  In  future  we  will  have  to  deal  with  cases 
in  which  similar  expansions  apparently  diverge,  we  will  not  tie 
up  the  development  of  our  method  of  the  construction  of  asymp¬ 
totic  approximations  with  any  proof  of  convergence. 

I  find  It  difficult  to  be  so  cavalier.  What  we're  supposed  to  do  is 
take  R 

R  =  -en2z2/2  *  -(en2a2/2)s1n21[sinwcosM+(l-e2)1/,2cosu)SlnM]2 


and  form  <R> 


If  a  term  of  <R>  depends  upon  an  element  that  has  a  secular  term,  and  the 
dependence  upon  this  element  is  period  (eg  the  sin  co, cos  <u  terms),  then  we 
are  to  say  that  these  are  long  periodic  terms.  If  the  dependence  is  not 
periodic,  then  these  are  secular  terms.  R-<R>  represents  the  short  periodic 
terms.  Does  <R>  have  a  secular  term?  Yes  if  you  say  rewrite  it  as 

<R>  .  -enVsir£i  (,.ecos2u) 


for  then 

Rsec  =  -en2a2sin2i/4 

2  2  2  2  2 
R£p  *  en  a  e  S1"n  icos 


No  if  you  don't.  Apparently  the  method  of  averaging  and  the  assignment  of 
terms  to  various  parts  of  R  depends  on  one's  mathematical  sophistication. 
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Ignoring  the  above  dilemnas  let's  turn  to  the  solution  for  the  orbital 

element  sets.  We're  supposed  to  obtain  the  short  period  changes  in  a,  e, 

etc.  by  using  R$p  in  Eqs.  (16).  This  is  fine  for  first  order  effects  except 

for  the  3R/3M  terms  because  the  operations  3/3M  and  / dt  don't  commute.  To 

evaluate  the  first  order  secular  terms  one  proceeds  in  a  similar  fashion. 

In  both  instances  one  regards  the  orbital  elements  appearing  on  the  right 

hand  sides  of  Eqs.  (16)  as  constants.  Again  this  is  fine  for  first  order, 

short  term  (times  <  P/e)  results.  Finally  one  is  supposed  to  obtain  the 

long  period  changes  by  an  exactly  analogous  procedure.  Since  w  varies  by 

2tt  on  a  time  scale  of  P/e,  the  periodicity  isn't  apparent  until  a  duration 
2 

of  order  P/e  »  and  regarding  a,e,i,etc.  to  be  constants  for  this  long  is 
of  dubious  validity. 

I  have  not  repeated  these  computations  here  because  they're  apparently 
not  rigorous,  they  don't  Illuminate  the  nature  of  the  motion  (short  term, 
secular, or  long  term)  and  they  can't  be  carried  out  in  an  unambiguous 
fashion  (if  I  can't  have  rigor,  I  at  least  demand  explicit  instructions). 

Had  I  produced  the  results  and  interpreted  them,  I  would've  drawn  some 
interesting,  and  false,  conclusions  regarding  the  motion  (remember  that 
this  is  an  exactly  soluble  problem  whose  solution  is  periodic).  Hence  I 
intend  to  turn  to  a  removal  of  part  of  the  problem  associated  with  perturba¬ 
tion  theory  in  celestial  mechanics--the  orbital  element  set  is  a  terrible 
basis  for  a  coordinate  system. 


C.  A  New  Representation 

The  above  "solution"  of  the  perturbation  equations  leaves  much  to  be 
desired— like  an  answer.  The  problem  Is  the  basis  of  the  representation 
a_-  (a,e,w,i ,  ft,T).  Go  back  to  the  original,  unperturbed  equations  of  motion 
(U  *  kr2/2), 

mr  =  F  =  -  VU  =  -kr 

The  general  solution  Is 

r  =  Acosnt  +  Bsinnt  ( 20) 

o 

where,  cf.  Eq.  (5),  n  =k/m.  If  we  add  a  perturbing  force  +ni7R  the  values  of 
A  and  B  will  no  longer  be  constant.  Let's  develop  a  perturbation  theory  for 
a  =  (A,jJ)  now.  Well  the  work  Is  all  done  and  Is  given  by  Eqs.  (14,15).  Once 
again  there  are  only  15  Independent  Lagrange  brackets  to  evaluate  since  21 
of  them  obviously  vanish.  These  are  of  three  types  too: 

Type  I:  both  a^.a^  are  one  of  the  elements  of  A  of  which 

3 

there  are  (2)  *  3  such, 

Type  II:  a^  is  one  of  the  elements  of  A. but  ak  is  one  of 
the  elements  of  B^  of  which  there  are  3*3=9  such,  and 
Type  III:  both  a^.a^  are  one  of  the  elements  of  B  of  which 
there  are  (2)  =  3. 

All  Type  I  and  III  brackets  vanish  as  do  the  "non- diagonal"  Type  II' s. 
The  only  non-zero  ones  are 

[AX,BX]  =  [Ay,By]  =  [AZ,BZ]  =  n 


30 


and  their  negatives.  The  simplicity  of  this  result  already  makes  one 
suspect  that  this  Is  a  very  nice  basis. 

The  perturbed  equations  of  motion  are 

mr  ■  -kr  -  (fe-k)z(0,0,l ) 

Let  fe  »  k(l+e)  with  |e|«l  as  before.  Then  R  s  -en2z2/2  or 
R  *  (Azcosnt  +  Bzs1nnt)2 

The  variation  of  parameters  equations  are  just 

•  •  •  • 

A  =  A  =  B  =  B  *0 
x  y  x  y 

Az  *  en(Azcosnt  +  Bzs1nnt)s1nnt  (21 ) 

Bz  *  -en(Azcosnt  +  Bzs1nnt)cosnt 

Much  simpler  It  would  be  difficult  to  Imagine.  Moreover  because  these 
equations  are  first  order,  linear,  ordinary  differential  equations  with 
polynomial  or  exponential  coefficients  Picard's  method  of  successive  approxi 
mat Ions  can  be  carried  out,  explicitly,  to  all  orders. 

Let  the  zero'th  order  approximation  be  denoted  by  Az(0),Bz(0).  These 
are  simply  related  to  the  Initial  values  for  z  and  z  In  the  unperturbed  case 

z (0)  =  Az(0),  z(0)  =  nBz(0) 

Using  primes  to  Indicate  the  order  In  Picard's  method  we  have,  successively 


Az  =  en[Az(0)cosnt  +  Bz(0)sinnt]sinnt 

Bz  =  -en[Az(0)cosnt  +  Bz(0)sinnt]cosnt 

Az  =  Az(0)  +  (e/2)Az(0)sin2nt  +  (e/2)Bz(0)(nt-sinntcosnt) 

Bz  =  Bz (0)  -  (e/2)Bz(0)sin2nt  -  (e/2)Az(0)(nt+sinntcosnt) 

^z  =  en(Azcosnt+Bzs1nnt)s1nnt 

Bz  =  -en(Azcosnt+Bzs1nnt)cosnt 

Az  =  Az  -  (e2/8)Az(0)[n2t2+(sinnt-2ntcosnt)sinnt]  + 

( e2/l 6 )BZ (0 ) ( 3si n2nt-4nt-2ntcos2nt ) 

Bz  ■  B^  +  (e2/8)Bz(0)[-n2t2+(3Slnnt-2ntcosnt)sinnt] 

+  (e2/16)Az(0)(sin2nt-2ntcos2nt) 

Az“=  n(A"zcosnt+Bzsinnt)s1nnt 
Bz" *  -  n(Azcosnt+Bzsinnt)cosnt 

etc.  We  see  short  period  terms,  secular  terms,  and  mixed  terms  in  the  above. 
Long  period  terms  will  never  appear  though. 

The  advantages  of  the  A,B  representation  over  the  orbital  element  set 
representation  should  be  clear  now.  Since  Picard's  method  can  be  rigorously 


applied  to  all  orders  all  we  need  to  do  is  sum  the  series  to  get  closed  form 


results  for  Az  and  Bz-  Of  course  for  this  problem  a  much  simpler 
available  to  us  for  there's  an  even  better  basis  than  the  A,f[  one. 
Define  a  and  8  via 

a  =  Azcosnt+Bzsinnt,  8  =  -Azsinnt+Bzcosnt 

They  satisfy 

a  =  n8»  8  =  -na(l+e) 


or 

a  =  -n2(l+e)a 

The  solution  to  this  can  be  written  down  by  inspection, 
a  =  Ccosvt  +  Dsinvt 

where  C  and  D  are  the  arbitrary  constants  of  integration  and 
v2  =  n2(l+e) 

Also  8  s  &/n  ■  (v/n)(-Csinvt  +  Dcosvt).  We  can  recover  Az  and  Bz 
Eqs.  (22), 

Az  =  acosnt-8sinnt,  Bz  =  asinnt+Bcosnt 
Or 


course  is 


(22) 


from 
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Az  =  C[cosntcosvt+(v/n)sinntsinvt] 

+  D[cosntsinvt-(v/n)s1nntcosvt] 

Bz  =  C[$inntcosvt-(v/n)cosntsinvt] 

+  D[sinntsinvt+(v/n)cosntcosvt] 

We  can  also  identify  C  and  D  as 

C  «  Az(0),  D  =  (n/v)Bz(0)  =  Bz(0)/(l+e)1/2 


(23) 


Complete  success--the  exact,  analytical  solution  of  the  perturbation  equations. 

It's  appropriate  now  to  reflect  on  our  accomplishments.  First  of  all  we 
found  a  basis  for  the  development  of  a  perturbation  theory  (a  set  of  arbitrary 
constants)  which  yield  a  much  simpler  set  of  variation  of  parameter  equations 
that  we  would've  expected.  Secondly  we  can  solve  this  by  actually  implement¬ 
ing  Picard's  method  of  successive  approximations.  Hence,  we  now  know  that  a 
solution  exists,  it's  continuous,  it's  unique,  etc.  Moreover,  for  this 
problem,  the  n'th  order  successive  approximation  in  Picard's  method  is  of 
order  e11.  This  coincidence  is  due  to  the  linearity  of  Eqs.  (21).  See  the 
next  section  for  a  counterexample  to  this  result  as  a  general  proposition. 

Thirdly  the  full  set  of  variation  of  parameters  equations  are  exactly  soluble 
in  closed  form.  This  reflects  the  fact  that  the  original  problem  was 
analytically  tractable  (since  a  =  2,  6  =  z/n)  and  the  choice  of  the  correct 
representation.  Finally  should  one  choose  to  expand  the  direct  solution  for 
Az,Bz  Eqs‘  (23)  in  a  Ta^or  series  in  e  one  will  recover,  order  by  order, 
the  successive  steps  of  the  Picard  scheme.  As  I  wrote  above,  a  complete  success. 
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V.  THE  ANHARMONIC,  ONE  DIMENSIONAL,  HARMONIC  OSCILLATOR 

In  this  Section  I  want  to  discuss  anharmonic  perturbations  and  since  the 
harmonic  oscillator  problem  separates  in  rectangular  coordinates,  one 
dimension  is  sufficient.  The  unperturbed  equations  of  motion  are 

2 

mx  =  -kx  ,  n  =  k/m 
whose  general  solution  is 

x  -  Acosnt+Bsinnt 

The  perturbed  equations  of  motion  are 

mx  =  -kx  -  nkx^ 

2  4 

so  R  =  -nn  x  /4  (I  could've  chosen  a  quadratic  form  for  the  anharmonic  term 
but  this  choice  simplifies  the  analysis  somewhat).  The  basis  for  the  pertur¬ 
bation  theory  is  a  =  (A,B).  The  only  non-zero  Lagrange  bracket  is  [A,B]  = 
-[B,A]  =  n.  The  variation  of  parameter  equations  are 

A  =  -(1/n)  8R/9B,  B  =  (l/n)9R/3A 
?  4 

Since  R  =  -nn  x  /4, 

R  =  (-nn2/4)  (Acosnt+Bsinnt)4 
or 

•  3 

A  =  nnsinnt  (Acosnt+Bsinnt) 

•  3 

B  =  -nncosnt  (Acosnt+Bsinnt) 

Let  A' ,B'  be  the  first  order  set  of  results  from  Picard's  method,  viz. 


A' (t)=A(0)+(n/32)[8AJ(0)(1-cos  7it)+3A^(0)B(0)(4nt-sin4nt) 
+24A(0)B2(0)sin4nt+B3(0)(12nt-8sin2nt+sin4nt)] 
B'(t)=B(0)-(n/32)[A3(0)O2nt+8sin2nt+sin4nt)+24A2(0)B(0)(l-cos4nt) 
+3A(0)B2(0)(4nt-sin4nt)+8B3(0)sin4nt] 

Clearly  the  next  approximation,  obtained  from 

•  o 
A"  =  nnsinnt(A'cosnt+B'sinnt) 

•  2 
B"  =  -nncosnt(A'cosnt+B'sinnt) 

results  in  powers  of  n  of  order  4  not  just  order  2.  Just  as  clearly  an 
anharmonic  perturbation  in  the  equation  of  motion  of  the  form  -nkxp~^  results 
in  A",B"  containing  terms  of  the  form  nP.  I  hope  that  this  result  clearly 
illustrates  the  difference  between  the  order  (in  the  sequence  of  successive 
approximations)  of  the  solution  and  the  functional  form  or  highest  power 
present  of  a  small  parameter. 

The  perturbation  equations  can  be  solved  analytically.  To  see  how 
define  a  and  B  via 

a  =  Acosnt+Bsinnt 
B  =  -Asinnt+Bcosnt 

Then,  because  of  the  condition  of  osculation. 


A  first  Integral  is 


a2  =  n2(L2-a2)  +  (nn2/2)(L4-ct4) 

if  the  initial  conditions  are  x(0)=L,  x(0)a0.  Other  initial  conditions  do 
not  yield  more  transparent  solutions  or  allow  me  to  make  additional  analytical 
points. 

The  solution  for  a  is 

a  =  Lcn[nt(l+nL2)1/,2.fe] 
fc2  *  L2/(2L2+2/n)  ,  fee[0,l] 

where  cn  is  the  cosine-amplitude  Jacobian  elliptic  function.  Since  3  =  a/n, 

6  =  -L(l+nL2)1/2sn[nt(l+nL2)1/2  ,  feldnCntO+nL2)17?  fe] 

where  sn  is  the  sine-amplitude  Jacobian  elliptic  function  and  dn  is  the 
delta-amplitude  Jacobian  elliptic  function.  The  solutions  for  A  and  B  are 
[m=nt(l+nL2)1/2] 

2  1/2 

A/L  =  cosntcn(m,fe)+(l+nL  )  sinntsn(m,fc)dn(m,fe) 

B/L  =  sinntcn(m,fe)-(l+nL2)^2cosntsn(m,fc),dn(m,fe) 

As  n-K) 
fe2^nL2/2 

2 

sn(m,fe)->sinm-(fe/2)  cosm(m-sinmcosm) 

2 

cn(m,fe)->cosm+(fe/2)  sinm(m-slnmcosm) 
dn(m,fe)-*-  1  -(fe2/2)2sin2m 


By  some  laborious  algebra  one  can  show  that  the  first  order  Taylor  series 

expansions  of  A  and  B  exactly  match  the  expressions  for  A'  and  B'[A(0)  = 

L,B  *0].  Clearly  the  statement  is  not  true  for  the  second  order  Taylor 

4 

series  since  A"  and  B"  contain  terms  proportional  to  n  •  I  don't  know 
if  the  fourth  order  Taylor  series  expansions  of  A  and  B  would  match  A"  and 
B" — I  have  only  a  finite  amount  of  patience  (I  would  bet  on  it  though). 

Lastly  let  me  mention  that  the  original  equations  of  motion  are  exactly 
soluble  too.  For  x(0)=L,  x(0)=0  the  solution  is 

x  =  Lcn(m,fe) 

-*■  Lcosnt  as  n-*0 

while  for  x(0)=0,  x(0)>0  it  is 

x  =  (L2+2/n)^2  sn(m,fe)/dn(m,fe) 

where  max(x)=L.  As  rr*  0  this  approaches  Lsinnt.  The  general  solution  is  not 
a  linear  combination  of  these  two  (it's  anharmonic)  nor  is  it  worth  displaying 


VI.  CONCLUDING  REMARKS 


Another  way  of  writing  the  general  solution  to  the  harmonic  oscillator 
(one  dimensional  for  the  moment)  is 

x  =  Ccos(nt+<J>) 

The  Lagrange  bracket  [C, <}>]=-  nC  and  the  perturbation  equations  are 
-  nC4>  =  3R/3C,  nCC  =  3R/3<j> 

But  R  is  time  dependent  for  either  anharmonic  or  anisotropic  perturbations. 
This  obviously  complicates  the  solution  and  I  haven't  explored  this  alterna¬ 
tive  in  depth.  Another  thing  I  haven't  done  is  derive  the  analog  of  Gauss's 
equations  for  non-conservative  perturbations.  As  an  example  if  there  is  a 
drag  term, 

mx  =  -kx  -  2myx 


then  the  general  solution  can  be  written  as  (I've  assumed  that  y  is  small 
2 

compared  to  n,  n  *  k/m  still) 

x  =  (Acosvt  +  Bsinvt)e”Yt 
2  2  2 

where  v  =  n  -y  .  It  would  be  interesting  to  pursue  this  problem  and  see 
where  else  the  usual  form  of  perturbation  theory  in  celestial  mechanics  let' 
us  down. 

An  even  more  interesting  generalization  would  be  the  simple  harmonic 
oscillator  subject  to  drag  and  an  external,  periodic  force.  The  equation  of 
motion  would  then  be  (say) 


mx  =  -kx  -  2myx  +  Fcos((ut+8) 


As  you'll  remember  the  steady  state  solution  Involves  a  resonance.  The 
phase  of  x  is  offset  from  that  of  the  external  force  by 

tani|f  =^4 
n  -a) 

and  the  amplitude  of  x  is  (F/m)/[(n2-uj2)2  +  4y2oj2]^2.  When  the  atmospheri 

2 

drag  is  small,  there  is  resonance  at  w  *  n-y  /n.  I  wonder  what  classical 
perturbation  theory  would  do  with  this. 
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